A numerical study of dynamics in thin hopper flow and granular jet
Wang Meng-Ke1, 2, Yang Guang-Hui1, 2, Zhang Sheng1, 2, Cai Han-Jie1, 2, Lin Ping1, 2, †, Chen Liang-Wen1, 2, Yang Lei1, 2
Institute of Modern Physics, Chinese Academy of Sciences (CAS), Lanzhou 730000, China
University of Chinese Academy of Sciences, Beijing 100049, China

 

† Corresponding author. E-mail: pinglin@impcas.ac.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11705256 and 11905272), the National Postdoctoral Program for Innovative Talents, China (Grant No. BX201700258), Young Scholar of CAS "Light of West China" Program for Guanghui Yang (Grant No. 2018-98), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA21010202).

Abstract

The dynamics of granular material discharging from a cuboid but thin hopper, including the hopper flow and granular jet, are investigated via discrete element method (DEM) simulations. The slot width is varied to study its influence on the flow. It is found the flow in the cuboid hopper has similarity with the flow in two-dimensional (2D) hopper. When the slot width is large, the flowrate is higher than the predicted value from Beverloo’s law and the velocity distribution is not Gaussian-like. For granular jet, there is a transition with varying slot width. For large slot width, there is a dense core in the jet and the variations of velocities and density are relatively small. Finally, the availability of continuum model is assessed and the results show that the performance with large slot width is better than that with small slot width.

PACS: ;81.05.Rm;
1. Introduction

There are several peculiarities distinguishing granular material from simple bulk solid, fluid or gas, such as internal friction, force chain, rotational degrees of freedom, and randomness in geometry.[13] Granular material can exhibit solid-like, fluid-like, and gas-like features under different boundary conditions and transitions often occur with a change of volume fraction.[4,5]

The granular materials discharging from hopper is an attractive issue for scientists and engineers. In the granular flow inside a cylindrical hopper, a famous rule is that the grains discharge constantly by gravity independent of filling height when the height is high enough (usually greater than the diameter of the hopper). In addition, the flow rate is independent of the diameter of hopper when it is larger than 2 times of the outlet diameter[6,7] and the outlet diameter is large enough to avoid any clogging.[8,9] This constant flowrate was observed by ancient people to invent the hourglass. The flowrate only depends on outlet and there is a scaling law,[10] which is studied in various cylindrical and square hoppers with conical or flat bottom. To explain this law, it is assumed to exist a free fall arch and the particles below the arch are falling freely.[11] The large-aspect-ratio hopper with a slot is rarely concerned in previous studies, which is proposed as a candidate configuration of high power target for neutron production.[12] Moreover, it is not very clear about the flow when the outlet size approaches to the hopper size.

As granular material is naturally discrete, the difference between granular and continuum fluid is one of most attractive issues for researchers. Because of the collisions, the motion of granular jet cannot be understood by analogy with free falling of one single particle. In Prado et al.’s work, it is interesting to find that the jet of powders will be thinning and this behavior is similar with incompressible fluid and is controlled by ratio of diameter of particle and diameter of funnel.[13] However, in simulation of impact on a target, the author suggested there is a significant difference between the granular jet and a fluid.[14] The experiments of the powder jet showed there are surface fluctuation of the jet, which can be explained by theory of capillary waves.[15] This instability will grow and clustering of powders will appear due to van der Waals interactions.[16,17] The clustering was also observed in simulation of cohesive granular jet.[18] Since the distinct difference between grains of micrometer level and powders, the behaviors of the jet of grains needs to be investigated deeply.

Recently, a phenomenological constitutive relation for plastic granular flow, which is called ‘μ (I) rheology’, was proposed. In Staron et al.’s work, the constant flowrate and pressure cavity above the outlet are successfully reproduced by employing μ (I) rheology to Navier–Stokes solver.[19] Besides usual hopper, they also obtain Beverloo scaling law in hopper with a lateral outlet.[20] This consistence indicates that granular flow in hopper can be treated as a Bingham fluid. It raises a question that if the continuum model is valid to simulate the granular jet.

This paper investigates the dynamics of spheres in a cuboid hopper with varying the slot size. The granular flow is divided into two parts: hopper flow and downstream granular jet. The profiles of velocity and volume fraction in these two parts are analyzed. The influence of outlet size on flow type of granular jet is studied. An incompressible fluid-like flow region with slight fluctuation is observed when outlet is large. The availability of continuum model in this flow is assessed.

2. Methods

The simplified model contains a cuboid hopper with a slot on the bottom. In the hopper there is a dense flow. There are many works about the cylindrical hopper but cuboid hopper is rarely concerned,[10,21] especially when the slot size approaches to the width of hopper. This hopper has a large aspect ratio (5:1) for statistical purpose. The spheres discharge from the slot and form a free-falling granular jet with gravity. In previous work for freely falling powders discharging from a container (< 100 μm), there is a clustering instability, boundary instability and turbulent behavior[15,17,22] and the spread angle of jet can be deduced.[23,24] But the constancy of the jet is still unclear. In this paper, the sphere diameter is 5 mm whose behavior is quite different with powders.[25] Due to the choose of size, in the simulations the effect of environment gas is neglected while it is obvious to powders.[26]

Due to the inhomogeneous and anisotropic properties, it is a big challenge to study granular matter theoretically and computationally,[2] and under some conditions the continuum approximation is not applicable anymore.[27,28] The discrete element method (DEM), closely related to molecular dynamics (MD), has become the most acceptable method for the study of granular flow since it was introduced by Cundall and Strack.[29] To simulate the large-scale granular flows, we developed the DEM code on GPUs,[30] which shows good validity and parallel efficiency.[3133]

Here we investigate the behavior of granular flow during the flat-bottomed hopper discharge through long narrow slot mounted on the bottom of the hopper. A system, consisting of 8 × 106 mono-disperse tungsten spheres (diameter d is 5 mm) is used and the material properties of the hopper is the same as spheres. The walls of the hopper are smooth and frictional. A Cartesian coordinate system is established and the center of the slot exit is set as the original point. The hopper dimensions in x and y directions are 200d and 40d respectively, resulting in an approximately 1000d height in z direction. The length of the narrow slot is fixed to be 200d and the slot width is a variable. The configuration of the flat-bottomed hopper is schematically shown in Fig. 1 and the material parameters are presented in Table 1. The time-step in DEM simulations is 5 × 10−7 s. Initially a random filling process is used to fill the hopper with the spheres until a random static packing state is reached (the volume fraction is about 0.58). Then the slot is opened and steady state flow is allowed to develop (generally 0.5 s after the flow starts) before acquiring the data used for subsequent analysis. In this paper the velocity and volume fraction fields of the hopper flow and granular jet are averaged in x direction. The derivations in the flow is also analyzed.

Fig. 1. Schematics of simplified model.
Table 1.

The parameters of the material properties in DEM simulations.

.
3. Results and discussion
3.1. Dynamics of hopper flow

In contrast with most fluids, granular materials display coexisting solid, liquid, and gaseous phases, which will produce a rich variety of complex flows. The hopper discharge through a bottom circular outlet has been extensively investigated both experimentally and computationally.[31,3335] To the best of our knowledge, the literatures on flow behavior through a slot of this type of flat-bottomed hopper, are very scarce and appear to be restricted to blockage problems.[36,37] Considering the shape of the flat-bottomed hopper, it can be inferred that spheres have no tendency to move in the x direction (as verified in Fig. 2). Thus, we just need to mainly study the flow pattern in yz cross section. In the following, we evaluate physical quantities from the height z = −100d to z = 40d. To obtain the spatial fields of these quantities, the yz cross section is divided into square cells of size 1.0d × 1.0d. Then in each cell the average/standard-deviation is performed over all the spheres of the center inside the cell and over all time points when the flow is steady. The profiles along the z direction are calculated as follows. In the z direction, the space is divided into bins of equal length Δz, generally set to 1.0d. For most quantities, averages are performed over these bins and over time points during the steady flow. But for the thickness of granular jet in the y direction, it is defined as the extension in that direction and the average is calculated over time.

Fig. 2. Spatial field plot of the time-averaged velocity in the x direction inside the hopper with the slot width D = 10d.

Firstly, we will analyze the simulation results of the flow inside the hopper. Granular hoppers have the well-known property of discharging their stored material at a constant rate, which depends on different parameters and is described by the Beverloo’s law.[10] To obtain the law of the flow rate in our case, we perform series of hopper discharges by varying the width of slot from D = 5d to D = 30d. Considering the observation that a large discrepancy between the measured flow rate and the value predicted by Beverloo’s law when the outlet size is significantly wider,[35] the slot width is limited to be less than 30d in our simulation. The flow rate is fitted by a modified Beverloo’s law for the hoppers with rectangular outlet: ,[38,39] where C0 and k are non-dimensional constants, ρf is the flowing bulk density at the outlet, g is the acceleration of gravity, d is the diameter of the spheres, L is the length of the outlet, and D is the width. The empirical term kd represents the “empty annulus” effect,[40] which indicates that no sphere center can approach within a distance of kd/2 to the outlet rim. Since L is fixed (L = 200d) and large enough in our case, LkdL and the formula can be rewritten as W = C(Dkd)3/2, where . The result shown in Fig. 3 indicates an agreement between the data and the fitted curve with parameters C = 22.47, k = 2.10. There is a discrepancy when D = 30d, which is in agreement with the experiments by Mankoc et al.[35] In a cylinder hopper, typical range for the value of k is usually between 1 and 2 and is independent of the sphere size.[41] In a two-dimensional (2D) case,[42] k = 2.17 for discrete simulations. The fitting results indicate that the law of flow rate in this configuration is similar to that in a 2D hopper.

Fig. 3. The flow rate as a function of the width (D) of the slot. The red line is a fit of Beverloo’s law. The inset shows a log–log plot.

Figure 4 shows the spatial vector field plot of the time-averaged vz in the yz plane. The vector plot shows that the velocity has a maximum at the center of the slot and the spheres gradually converge to the slot. Two small regions in the left and right corners made by the sidewalls and bottom wall remain stagnant.[43] As can be seen in Fig. 4, the flow in the upper hopper is a plug-flow. In addition, we observe the nonzero slip velocity at sidewalls, which depends on the friction properties of the contacts between spheres and walls.[44,45] Figure 5 shows the spatial field plot of standard deviation of the velocity magnitude in the yz plane. It shows a strong variation near the slot rim, which is perhaps attributed to frequent collisions between spheres and the slot rim. This phenomenon agrees with the fact that grains close to solid boundaries in granular flows are hotter in terms of granular temperature.[4648] In this paper, we use a Voronoi tessellation to calculate the volume fraction on a local scale, which uniquely assigns a polygonal volume to each sphere. The volume fraction is defined as the ratio of the sphere volume to the polygonal volume, and is calculated using the algorithms proposed by Rycroft.[49]

Fig. 4. Spatial vector field plot of time-averaged vz in the yz plane inside the hopper with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d. The arrow heads indicate the velocity direction. The length of arrows and the color scale indicate the velocity magnitude.
Fig. 5. Spatial field plot of standard deviation of vz in the yz plane inside the hopper with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = v30d.

Figure 6 shows the spatial field of time-averaged volume fraction. It shows that there is volume fraction with a constant value of 0.55 in plug-flow region, while there is a sharp transition to low volume fraction (0.45–0.35) near the slot, due to the shear dilation.

Fig. 6. Spatial field plot of the time-averaged volume fraction inside the hopper with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d.

The downward velocity profiles inside the hopper with D = 10d at different heights are plotted in Fig. 7. The profiles near the slot have a parabolic shape. According to the rheology of granular flow, this region is classified as the inertial (rapid flow) regime. When z > 10d, the profiles indicate that the flow is a plug flow and is located at the quasi-static regime. In Refs. [50,51] a continuum kinematic model containing only one unknown parameter is proposed, and it can provide an accurate prediction for the description of the velocity profiles in a quasi-2D flat-bottomed bottom hopper with a centric discharge. Based on the following two assumptions: the constitutive law and the boundary condition, this model gives a Gaussian-like profile for downward velocity

where W is the volumetric flow rate per unit thickness of the hopper and b is referred to as the ‘diffusion length’ by Choi et al.[52,53] Then we compare the velocity profiles with the prediction of this model by using similar methods in Refs. [52,54] where the non-linear least square is used to find the values of and b. In Fig. 8, we show the downward velocity profiles for D = 5d, 10d, 20d, 30d at two heights z1 = 2d and z2 = 5d along with the fit to the kinematic model. For small slot width (D = 5d). The profiles at z1 and z2 are best fitted with b = 1.2d and b = 2.0d respectively. For intermediate slot width (D = 10d), the diffusion lengths b = 2.8d and b = 3.0d are the best fit for the profiles at z1 and z2. The best fit lines for small and intermediate slot width show good agreement with velocity profiles. In addition, the diffusion length b becomes larger when z increases, in accordance with some previous reports.[5255] However, for large slot width (D = 20d and 30d), the shape of velocity profile at z1 is obviously deviated from the Gaussian distribution as it has a flat sharp at the center. And for large slot width, the best fitting values of b (15.0d and 8.7d for z1 and z2 respectively) are out of the normal range (< 4d).[50,55,56]

Fig. 7. The time-averaged profile of vz as a function of y coordination inside the hopper with D = 10d at different heights. The vertical black dashed line is drawn for visual clarity.
Fig. 8. The time-averaged profile of vz for D = 5d, 10d, 20d, 30d at two heights z1 = 2d and z2 = 5d. The fit lines using kinematic model for heights z1 = 2d and z2 = 5d are also shown. The vertical black dashed line is drawn for visual clarity.
3.2. Dynamics of granular jet

Granular jets falling under gravity exhibit astonishingly liquid-like appearance.[57] In recent papers,[15,16,58,59] such jets were investigated, where the main interest is focused on describing the surface properties of these jets like the clustering instability due to cohesive force. Similar to liquid jets, these jets became progressively thin as they flow downward, with a roughly constant density implying incompressible-like feature. However, Prado et al.[47] observed a clear transition from compressible to incompressible granular jets in their experiments, which depends solely on the aspect ratio between the diameter of spheres and the diameter of the outlet. In the following section, we will investigate the flow behaviors of the free-falling granular jet emanating from the long slot. It is interesting that, we find a similar “compressible-incompressible” transition in our simulations by varying the slot width D.

Figure 9 shows snapshots of granular jets falling under gravity for different slot width D. For small slot width (D = 5d, 7d, 10d), the granular jets spread out when they propagate downward. However, for intermediate slot width (D = 15d) the jet spread out slightly, and for large slot width (D = 20d, 30d) the jet thins at first and then scarcely spreads out. Since the displacement and velocity of each sphere is known individually in the DEM simulations, many quantities are easy to obtain. The thickness (L) of granular jet, which is defined as a time-averaged position of outmost particles, is measured to characterize the jet dispersion. To describe the contraction or expansion of this jet, figure 10 illustrates the rescaled thickness (L/D) of jets as function of z for different slot width D. As we can see, the degree of granular jet dispersion decreases with the slot width. For small slot width (e.g., D = 5d), the granular jet disperses strongly in the region far from the slot, and even at z = −100d the thickness (L) is about 6 times of the slot width (D). An interesting phenomenon is that for large slot width (D = 20d) the thickness of jet decreases slightly at first and then increases gradually along the z axis. The above phenomenon are basically in accordance with the snapshots shown in Fig. 9. In addition, the thickness of granular jet as a function of time at different heights for D = 10d is illustrated in Fig. 11. It shows that the rescaled thickness (L/D) fluctuates with time but shows a small variation, which indicates that the flow rate of the jet flow is rather steady.

Fig. 9. Snapshots (in yz cross section) of granular jets near the slot of the hopper for different slot widths D (from left to right 5d, 10d, 15d, 20d, 30d).
Fig. 10. The rescaled thickness of granular jet as a function of z for different slot widths D. The vertical black dashed line indicates no dispersion.
Fig. 11. The rescaled thickness of granular jet as a function of time at different heights for slot width D = 10d.

As mentioned above, compared with physical experiments, DEM simulations have the ability to acquire the information within the jets. Figure 12 shows the spatial vector field plot of the time-averaged vz in the yz plane of granular jets. When D = 30d, there is a valley of vz in horizontal direction, which is consistent with the velocity profile inside the hopper (see Fig. 8). Compared with Fig. 9, there is no contraction at z = −10d which is due to the infrequent scattered particles outside the jet. The spatial vector field of standard deviation of the velocity magnitude in the yz plane is displayed in Fig. 13. It can be seen that the relatively high variations appear to be near the edge of the slots, which is similar to the findings shown in Fig. 4.

Fig. 12. Spatial vector field plot of the time-averaged vz in the yz plane of granular jets with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d. The arrow heads indicate the velocity direction. The length of arrows and the color scale indicate the velocity magnitude.
Fig. 13. Spatial field plot of standard deviation of vz in the yz plane of granular jets with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d.

The center-line downward velocity at different heights is also studied and the profile for D = 10d is presented in Fig. 14. If there is no collision between the spheres at the center-line, wherein the sphere is accelerated by gravity and the velocity grows with z as

It can be seen that the sphere velocity has nearly no difference with the velocity predicted by the free-falling model, which indicates there are little collisions between the spheres during the falling phase of granular jets.

Fig. 14. The vz and volume fraction as a function of z for the granular jet with D = 10d. The solid line shows the velocity of free-falling single sphere.

The time-averaged profiles of downward velocity at different heights for D = 10d are displayed in Fig. 15. As is expected, the velocity increases with the growth of the distance from the slot exit. And the spheres at the center region have larger velocity than those at the boundary. As reported in Refs. [60,61] the Gaussian profile can be used to describe to the far-field downward velocity of the granular jet. However, a recent experiment[62] suggested that the downward velocity distribution at the exit is a parabolic function, and can be written as with considering the vertical free fall motion. However, the differences between the fit lines in Fig. 15 and the data at the boundary are attributed to the fact that the horizontal velocities of these spheres, leading to the expanding, are nonzero.

Fig. 15. The profile of vz in granular jet with D = 10d at different heights. The fit lines using parabolic model are also shown. The vertical black dashed line is drawn for visual clarity.

Figure 16 shows the spatial field plot of the time-averaged volume fraction of granular jets. As we can see, the volume fraction at the center is higher than that at the boundary. For small and intermediate slot width, the volume fraction at the center decreases significantly along the z-axis due to the dispersed feature of the jets. However, for large slot width, there is nearly no dispersion and the velocity of sphere increases along the z axis, causing slight decreases in the volume fraction. The volume fraction of granular jet as a function of z for different slot width D is presented in Fig. 17. A notable decrease of the volume fraction along the z axis is observed for small slot widths, while a relatively slow decrease is observed for the high slot widths. These observations are in accordance with the field plots shown in Fig. 16. The spatial field plot of the standard deviation for the volume fraction is displayed in Fig. 18. It is worthwhile to note that the large variations appear to be located near the periphery of the slot, and spread downward along the z axis.

Fig. 16. Spatial field plot of the time-averaged volume fraction of granular jets with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d.
Fig. 17. The time-averaged volume fraction of granular jet as a function of z for different slot widths D.
Fig. 18. Spatial field plot of standard deviation of the volume fraction in the yz plane of granular jets with the slot width (a) D = 5d, (b) D = 10d, (c) D = 20d, (d) D = 30d.
3.3. Availability of continuum method

As mentioned above, the continuum model was successfully to reproduce some important features of hopper fluid. But a quantificational evaluation must be performed for the hopper flow and to our knowledge, this model has not been used to describe the granular jet.

We use ANSYS Fluent software[63] to simulate the continuum counterpart of the granular jets. Following Staron et al.,[19] we set = min (μ P/D2ηmax, where ηmax = 1000 Pa⋅s and μ is the effective coefficient of friction of the granular flow: μ = μs + Δ μ/(1+I0/I), where μs = 0.32, Δ μ = 0.28, and I0 = 0.4 are constants. P is the local pressure and D2 is the second invariant of the strain rate tensor D: is a dimensionless parameter: , where ρ is the density of spheres. The simulation is simplified to a 2D model, and viscous model is set to standard K-Epsilon. A no-slip boundary condition is imposed at the wall near which the standard wall functions is used.

We simulate the hopper flow with D = 10d and D = 30d. The flowrates are consistent with the results of DEM simulation. The spatial profiles of vertical velocity of hopper flow and granular jet are plotted in Fig. 19. Because the constant volume fraction is used in the continuum model, the jet becomes thinner even when D = 10d as can be seen in Fig. 19. Furthermore, we plot the spatial profile of ratio of error, which is define as , where vz is the local vertical velocity in DEM simulation and is from simulation of continuum model. Figure 20 shows there is an obvious discrepancy in velocity field with small outlet (D = 10d) but for large outlet (D = 30d) the performance of continuum model is quite better. In granular jet with D = 30d, except the boundary, the discrepancy is no more than 20%.

Fig. 19. Spatial profiles of vz in hopper flow with (a) D = 10d, (b) D = 30d. Spatial profiles of vz in granular jet with (c) D = 10d, (d) D = 30d.
Fig. 20. Spatial profiles of ratio of error in hopper flow with (a) D = 10d, (b) D = 30d. Spatial profiles of ratio of error in granular jet with (c) D = 10d, (d) D = 30d.
4. Conclusion

In this paper, we simulate the spheres discharged from a cuboid, flat-bottomed hopper with different slot widths. The granular flow inside the hopper is investigated by plotting the spatial profiles of velocity and volume fraction. It shows when the slot width is small, the velocity is nearly constant with the time going by and the distribution can be predicted with a Gaussian-like function. However, the prediction fails when the slot is wide and the flowrate is higher than the predicted value from Beverloo’s law. There is a high fluctuation of vertical velocity at the rim of outlet. The volume fraction in the hopper flow is also constant but decreases above the outlet. In the granular jet beneath the hopper, the volume fraction will drop with height rapidly because of the dispersion of granular jet if the slot width is small. In contrast, there is a dense core in granular jet if the slot width is large. The fluctuations of velocity and volume fraction in this core are relatively small. The recent continuum model is also investigated and the results indicate a potential availability.

The dynamics of granular jet can be described by a free-falling function, which indicates a low energy dissipation. Further work needs to be done to confirm this point. Another question is that what will happen if slot width approaches the hopper width? In this paper it shows the invariable volume fraction limits the applicability of continuum model, which should be improved in the next step.

Reference
[1] Jaeger H M Nagel S R 1992 Science 255 1523
[2] de Gennes P G 1999 Rev. Mod. Phys. 71 S374
[3] Chuang Zhao C B L Bao Lin 2018 Chin. Phys. 27 104501
[4] Brodu N Delannay R Valance A Richard P 2015 J. Fluid Mech. 769 218
[5] Wang W G Zhou Z G Zong J Hou M Y 2017 Chin. Phys. 26 044501
[6] Nedderman R M Tuzun U Savage S B Houlsby G T 1982 Chem. Eng. Sci. 37 1597
[7] Zhang S Lin P Yang G Wan J F Tian Y Yang L 2019 Chin. Phys. 28 018101
[8] To K Lai P Y Pak H K 2001 Phys. Rev. Lett. 86 71
[9] Ma L D Yang G H Zhang S Lin P Tian Y Yang L 2018 Acta Phys. Sin. 67 044501 in Chinese
[10] Beverloo W A Leniger H A Vandevelde J 1961 Chem. Eng. Sci. 15 260
[11] Brown R L 1961 Nature 191 458
[12] Cai H J Yang G Vassilopoulos N Zhang S Fu F Yuan Y Yang L 2017 Phys. Rev. Accel. Beams 20 023401
[13] Prado G Amarouchene Y Kellay H 2013 Europhys. Lett. 102 198001
[14] Sano T G Hayakawa H 2012 Phys. Rev. E 86 1
[15] Amarouchene Y Boudet J F Kellay H 2008 Phys. Rev. Lett. 100 4286
[16] Royer J R Evans D J Oyarte L Guo Q Kapit E Mobius M E Waitukaitis S R Jaeger H M 2009 Nature 459 1110
[17] Mobius M E 2006 Phys. Rev. 74 051304
[18] Ulrich S Zippelius A 2012 Phys. Rev. Lett. 109 166001
[19] Staron L Lagree P Y Popinet S 2012 Phys. Fluids 24 47
[20] Zhou Y Lagree P Y Popinet S Ruyer P Aussillous P 2017 J. Fluid Mech. 829 459
[21] Peng L Xu J Zhu Q S Li H Z Ge W Chen F G Ren X X 2016 Powder Technol. 304 218
[22] Uchiyama T Naruse M 2006 Chem. Eng. Sci. 61 1913
[23] Cao W G Liu H F Li W F Xu J L 2014 Fuel 115 17
[24] Prado G Amarouchene Y Kellay H 2011 Phys. Rev. Lett. 106 198001
[25] Caretta O Loveridge P O’Dell J Davenne T Fitton M Atherton A Densham C Charitonidis N Efthymiopoulos I Fabich A Guinchard M Lacny L J Lindstrom B 2018 Phys. Rev. Accel. Beams 21 101005
[26] Ma L Zhang X Zhang S Lin P Zhang Y Liu W Sun J Zhu Y Xiao R Yang G Tian Y Yang L 2018 Nucl. Eng. Des. 330 289
[27] Schaeffer D G 1992 Proc. Roy Soc. Lond. Mat. 436 217
[28] Ostoja-Starzewski M 1993 Mech. Mater. 16 55
[29] Cundall P A Strack O D L 1979 Géotechnique 29 47
[30] Tian Y Zhang S Lin P Yang Q Yang G Yang L 2017 Comput. Chem. Eng. 104 231
[31] Zhang S Lin P Wang C L Tian Y Wan J F Yang L 2014 Granul. Matter 16 857
[32] Lin P Zhang S Qi J Xing Y M Yang L 2015 Physica A: Statistical Mechanics and its Applications 417 29
[33] Tian Y Lin P Zhang S Wang C L Wan J F Yang L 2015 Adv. Powder Technol. 26 1191
[34] Beverloo W A Leniger H A van de Velde J 1961 Chem. Eng. Sci. 15 260
[35] Mankoc C Janda A Arévalo R Pastor J M Zuriguel I Garcimartín A Maza D 2007 Granul. Matter 9 407
[36] Langmaid R N April H E 1957 J. Inst. Fuel 166
[37] Davies C E Desai M 2008 Powder Technol. 183 436
[38] Myers M Sellers M 1978 University of Cambridge, Cambridge
[39] Anand A Curtis J S Wassgren C R Hancock B C Ketterhagen W R 2008 Chem. Eng. Sci. 63 5821
[40] Brown R L Richards J C 1970 Principles of powder mechanics: essays on the packing and flow of powders and bulk solids Pergamon Press 495
[41] Hirshfeld D Rapaport D C Eur. Phys. J. E 4 193
[42] Staron L Lagrée P Y Popinet S 2014 Eur. Phys. J. E 37 51
[43] Nedderman R M 1995 Chem. Eng. Sci. 50 959
[44] Artoni R Santomaso A C Go’ M Canu P 2012 Phys. Rev. Lett. 108 238002
[45] Artoni R Richard P 2015 Phys. Rev. Lett. 115 158001
[46] Andreotti B Yoël F Olivier P 2011 Les milieux granulaires EDP Sciences
[47] Prado G Amarouchene Y Kellay H 2013 Europhys. Lett. 102 24006
[48] Sano T G Hayakawa H 2012 Phys. Rev. 86 041308
[49] Rycroft C H Grest G S Landry J W Bazant M Z 2006 Phys. Rev. 74 021306
[50] Nedderman R M Tüzün U 1979 Powder Technol. 22 243
[51] Tüzün U Nedderman R M 1979 Powder Technol. 24 257
[52] Balevicius R Kacianauskas R Mroz Z Sielamowicz I 2011 Adv. Powder Technol. 22 226
[53] Jaehyuk C Arshad K Martin Z B 2005 J. Phys.: Condens. Matter 17 S2533
[54] Balevicius R Sielamowicz I Mroz Z Kacianauskas R 2011 Powder Technol. 214 322
[55] Medina A Córdova J A Luna E Treviño C 1998 Phys. Lett. 250 111
[56] Mullins W W 1979 Powder Technol. 23 115
[57] Jens E Emmanuel V 2008 Rep. Prog. Phys. 71 036601
[58] Möbius M E 2006 Phys. Rev. 74 051304
[59] Ulrich S Zippelius A 2012 Phys. Rev. Lett. 109 166001
[60] Ansart R Ryck A d Dodds J A 2009 Chem. Eng. J. 152 415
[61] Ansart R de Ryck A Dodds J A Roudet M Fabre D Charru F 2009 Powder Technol. 190 274
[62] Janda A Zuriguel I Maza D 2012 Phys. Rev. Lett. 108 248001
[63] Fluent I 2012 Fluent Inc. Canonsburg PA, USA